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Abstract 


Introduction: Both in-vitro and in-vivo observations and experimental research have shown that 
the radial distribution of red blood cells in blood flowing in arteries is, under certain 
circumstances, not uniform. Straight segments of small to medium caliber arteries show an axial 
accumulation of blood cells. A field of apparent variable viscosity appears because of this 
distribution pattern. Recent research in computational bio-fluid mechanics shows that fluid 
velocity profiles and wall shear stress in small arteries are affected when the macro-rheological 
Newtonian and Non-Newtonian blood viscosity mathematical models are compared. Commonly 
used non-Newtonian models don't consider non-uniform radial distribution of blood cells, and 
consequently the spatial variation of viscosity. Moreover, variations in shear rates modify fluid- 
to-wall mass transport processes in arteries, related with the physiology and pathophysiology of 
arteries. 


Objectives: To consider this spatially variable viscosity as it is related to and affects 
hemodynamics and fluid-to-wall mass transport processes in arteries, emphasizing an 
approximate analytical approach contrasted with digital simulations. 


Materials & methods: One phase incompressible and isotropic theory of viscous suspensions of 
non-Brownian particles is used to describe blood flow in in straight segments of small arteries. 
A mathematical model of a spatially variable hematocrit is constructed from given axisymmetric 
velocity patterns in straight segments of small arteries under quasi-stationary conditions. 
Velocity profiles, pressures, local viscosity fields, shear rates and wall shear stresses are then 
obtained using finite elements digital simulation, in straight segments of small arteries under 
pulsated flow conditions. The effects of the resultant hemodynamics on the parameters of fluid- 
to-wall mass transport processes are estimated from wall shear rates by means of diffusive 
boundary layer theory. 


Results: In the considered global (tube) haematocrit range (35-55%), obtained shear rate profiles 
and transport parameters in fluid adjacent to the wall are affected by the macro-rheological 
model used and by the local haematocrit field according to the micro-rheological model 
proposed. 


Conclusions: Considering spatial variations of viscosity in bio-fluid computational mechanics 
simulations, may modify significantly local hemodynamics and fluid-to-wall transport processes 
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as well as pumping power needed to maintain pulsated flow in a small artery, both with 
Newtonian and Non-Newtonian rheology. 


Keywords: computational bio-fluid mechanics, rheological models of blood, models of axial 
accumulation of blood cells, mass transport processes, diffusion boundary layer. 


(A) Introduction 


Already in the 18th century, Haller and Spallanzani noticed during careful in-vivo 
observations that blood cells tended to accumulate towards the axis of the vessel. In 
1941, Miiller observed in-vitro the same behavior, working with a straight capillar tube 
and a laminar steady flow of a glycerol solution with (same density) suspended micro 
rubber discs (Taylor 1955; Bayliss 1959). Since 18th century, both in-vitro and in-vivo 
observations and experimental research have shown that the radial distribution of red 
blood cells in blood flowing in arteries is often not uniform (Rampling 2007). Besides 
the always present almost cell-free and slender plasma layer adjacent to the 
endothelium, there is a variation of the distribution of blood cells inside the lumen of 
the vessel. The accumulation pattern, if present, depends on the pattern of blood flow. 


At first, blood viscosity measurements in cylindrical tube viscometers (with capillary 
tubes) were done applying Poiseuille formula for a homogeneous Newtonian fluid: 


_ am-R*- AP 
8-L-O 


Here R is tube's radius, L is tube's length, Q is blood flow, and AP is pressure 


difference between tube ends. Tube length L is big enough, so entrance length effects 
can be neglected. 
The viscosity 4 should be a blood parameter independent of tube geometry. However, 


experimental measurements show that there is a relation between apparent viscosity and 
tube radius. This is the so-called Fahraeus-Lindqvist effect: the apparent viscosity 
decreases when the tube radius decreases (Ross-Ethier and Simmons 2007; Cokelet and 
Meiselman 2007a). This gradual decrease in viscosity with decreasing tube size 
becomes detectable with tubes less than 500um in internal radius. This reduction 
becomes significant with tubes less than 200um in internal radius an up to 50 um. 


There are two classical macro rheological models that explain Fahraeus-Lindqvist 
effect: wall effect model and sigma effect model. 

Wall effect model splits the Newtonian fluid in two phases, a cell free marginal plasma 
layer, adjacent to the inner wall surface, and a cylinder of uniformly packed red blood 
cells. Apparent viscosity is given as a function of the quotients between plasma layer 
thickness (with a thickness of around 3 um) and tube radius and between plasma 
viscosity and the viscosity of uniformly packed red blood cell cylinder (Cooney, 1976). 
Sigma effect model splits the Newtonian fluid in a finite number of concentric lamina, 
related to the discrete nature of red blood cells. Each lamina is completely un-sheared 
internally, and the velocities of adjacent lamina are different. Apparent viscosity is 
given as a function of plasma viscosity and the quotient between lamina thickness and a 
uniform whole blood viscosity (Cooney, 1976). Wall effect and sigma effect models are 
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usually combined to give a more complete formula for the apparent viscosity (Cooney, 
1976, Sudrez-Antola, 1980). 

Besides, several empirical correlations between blood viscosity, tube radius and global 
hematocrit were proposed (Baskurt et al. 2007). 


Experimental evidence (Pearl, 1961; Cooney, 1976; Baskurt and Meiselman 2007) 
shows that straight segments of small to medium (from several hundred of um to mm) 
caliber arteries present an axial accumulation of blood cells. At low enough shear rates, 
some of them travel aggregated in rouleaux and others travel free. A field of apparent 
variable viscosity appears because of this distribution pattern, growing towards vessel 
axis. To consider this spatial variation, and to make explicit its effect on apparent 
viscosity, the abovementioned combination of the classical wall effect and sigma effect 
models is not enough: a different mathematical model must be introduced, even if the 
fluid is considered as Newtonian. However, below shear rates of about 50 s™! apparent 
blood viscosity decreases with increasing shear rates (blood thinning), so non- 
Newtonian rheological models must be used in this case (Cokelet and Meiselman 
2007b). Experimental evidence related to the effects of blood cell migration on local 
viscosities, associated with complex blood flows (related with vessel curvature, 
ramifications, bifurcations, and changes in cross sections) found in larger vessels (from 
mm to cm) is fairly scarce. The elucidation of these effects is important because blood 
viscosity is a key parameter that modulates hemodynamic forces and fluid-to-wall mass 
transport processes. 


One of the main consequences of this modulation is an effect of mechanical 
transduction at the endothelial level. Endothelial cells that line the inner surface of the 
artery sense and responds to shear stress. Shear stress patterns affect endothelium 
permeability and the release of signaling molecules (Nichols et al. 2011). The 
magnitude and pattern of shear stress are related with the physiology and 
pathophysiology of arteries (Nichols et al. 2011), and are characterized by several 
indexes like OSI, WSSG, etc. (Chandran et al. 2012). Moreover, hemodynamics 
influence fluid-to-wall mass transport processes in arteries and these processes are 
related also with arterial physiology and pathophysiology, such as atherosclerotic 
plaque location and development (Ross-Ethier 2002; Tarbell and Qiu 2006; Coppola 
and Caro 2009; Caro et al. 2012). 


Different substances like oxygen and low-density cholesterol (LDL) play a key role in 
the development and progression of atherosclerotic disease. These substances travel 
from blood to the artery wall. The transport of some of these substances, like oxygen 
and other small molecular weight substances may be limited by the convective- 
dispersive flow from bulk blood to endothelial surface. Between certain bounds, this 
flow increases with the local fluid shear rate adjacent to the wall. The transport of LDL 
and other high molecular weight substances seems to be limited mainly by the 
endothelium. As consequence, as already mentioned above, their flow is modulated by 
shear stress at the fluid-structure interface. 


If blood is modeled as a homogeneous incompressible Newtonian fluid, then shear 
stress is always proportional to shear rate under any flow conditions. But if flow is 
modeled as a homogeneous incompressible non-Newtonian fluid, in general the 
proportionality between shear stress and shear rate vanishes and as consequence 
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mechanical transduction processes doesn’t parallel fluid-to-wall mass _ transport 
processes. Recent research in computational bio-fluid mechanics shows that fluid 
velocity profiles and wall shear stress in small arteries are significantly affected when 
the Newtonian and the non-Newtonian mathematical models of blood viscosity (for 
homogeneous fluids without axial accumulation of blood elements or cell free plasma 
layer) are compared. Kanaris et al. (2012) compared the results of CFD simulations 
obtained applying a Casson rheological model with the corresponding results obtained 
with a Newtonian model for a vessel bifurcation. In both cases they considered the 
dependence of viscosity parameters on the global hematocrit. 


If there is a spatially and temporally variable hematocrit field, even in the Newtonian 
case, without further analysis, the proportionality between shear stress and shear rate 
should not be taken as granted. 

By construction, neither the classical Newtonian nor the commonly used non- 
Newtonian phenomenological models (like Quemada, Carreau, Bingham, Casson, and 
Potential Law), consider a possible non-uniform radial distribution of blood cells, and 
consequently the corresponding spatial variation of viscosity due to particle migration in 
the flow field (Woodcock 1976, Yilmaz and Gundogdu 2008). 

However, to assess the effects of this variable viscosity as it is related to and affects 
hemodynamics and mass transport processes, it is necessary to consider spatial variation 
of parameters (mainly due to the migration of blood cells) in blood rheological models 
used in computational bio-fluid mechanics. 


Particle migration is a common phenomenon in flowing suspensions of both, rigid and 
deformable particles (Davis 1996, Shauly et al. 2000, Scott, 2005, Zhao et al. 2012, van 
Dinther et al. 2013). In general, migration occurs from regions of high shear rates to 
regions of low shear rates'. This fact can be explained, from a non-equilibrium 
thermodynamics point of view, applying the principle of minimum entropy production 
(Quemada 1977, Sudrez-Antola 1980, Glaser 2001, Dillon 2012). While the 
thermodynamic approach allows an interpretation of particle migration using a common 
framework, the mechanisms that effectively produce the migration of the suspended 
particles differ depending on the concentration of the suspension, the size, the shape and 
the degree of deformability of the suspended particles, buoyancy effects, and flow 
regime, amongst other properties (Sudrez-Antola 1980). 


The main migration mechanism in the case of concentrated suspensions of deformable 
non-Brownian particles, such as red blood cells, is the so-called hydrodynamic diffusion 
(Phillips et al. 1992, Davis, 1996, Biasetti et al. 2014). The formula for the suspended 
particles flux density, derived by Pillips et al. in their phenomenological model of 
hydrodynamic diffusion, can be re-written (van Dinter et al. 2013) as the product of a 
positive scalar conductance field by minus the gradient of another scalar field, known as 
the migration potential. As consequence of this result, suspended non-Brownian 
particles always move towards regions of lower migration potential. Equilibrium is 


1 This remark applies to a dominant population, in the sense that the interaction between its 
particles and between the carrier fluid and the particles determine the main characteristics of the 
flow and in particular its energy dissipation. This is the case of the population of red cells in 
blood. However, the behavior of minority populations, like platelets, white blood cells and 
eventually drug-delivering particles, may be different and even opposite: instead of migrating 
towards the conduit axis, the move towards the wall (Baskurt et al. 2007, Zhao et al. 2012). 
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attained when the migration potential is the same in all points of the fluid. 

The conductance field is proportional to the square of the local particle concentration 
(the local hematocrit if the fluid is blood), the square of a characteristic particle size, 
and a positive measure of the local shear rate in the fluid. The migration potential field 
is the logarithm of the product of the local particle concentration, the above-mentioned 
measure of the local shear rate and a power of the local fluid viscosity. 


Recently, Gidaspow and Huang, with the intention of attaining a comprehensive micro 
rheological theory of Fahraeus-Lindqvist effect, developed a kinetic theory of shear 
induced particle migration in the framework of a two-phase fluid model for blood flow 
in narrow straight cylindrical tubes (Gidaspow and Huang 2009). 

In this case the migration of the non-Brownian particles towards the axis of the flow is 
driven by a gradient in a local granular temperature that may be defined for a population 
of suspended particles in a flow. This temperature decreases radially from the wall to 
the axis. 

Defined as the rate of momentum transfer per unit area due to the interactions between 
the deforming, translating, and rotating suspended particles, a local collision pressure 
(roughly proportional to the above-mentioned granular temperature) may be introduced 
(Sudrez-Antola 1980). This collision pressure (closely related with the migration 
potential of the hydrodynamic diffusion) appears in a Fokker-Planck equation that 
approximately describes the time and spatial evolution of the local hematocrit. The 
spatially variable viscosity of the suspension can be calculated by statistical mechanics 
methods (Gidaspow and Huang 2009). 

As far as the authors know, this kinetic approach has not been applied yet to study 
migration of deformable particles in complex flows. 


Biasetti, Spazzini, Hedin and Gasser (Biasetti et al. 2014) published the results of a 
CFD research that considers shear induced migration of blood cells in complex flow 
fields like the ones found in medium to large arteries, working with a non-Newtonian 
rheological model with realistic values of the global hematocrit. They combined, in 
single phase non homogenous CFD model, Phillips phenomenological model of 
hydrodynamic diffusion with Quemada’s non-Newtonian phenomenological viscosity 
model. With the obtained CFD model, an in-silico simulation of the behavior of the 
formed elements was done: in a normal abdominal aorta, a normal carotid bifurcation, a 
stenotic carotid bifurcation and an abdominal aneurysm. The results of the digital 
simulations suggest a significant blood cell migration in the first three cases, but not in 
the abdominal aneurysm. 


(B) Objectives 


The general purpose of the present work is to study some aspects of the influence of a 
spatially variable viscosity on hemodynamics and fluid-to-wall mass transport processes 
in arteries, using a mixed analytical and numerical framework and emphasizing the 
analytical approach. To ease the application of analytical methods, a simplified model 
of axisymmetric flow in a straight segment of artery is considered under quasi- 
stationary conditions, so that the shear stress in the lumen is known and independent of 
the rheological model. 


Under these restrictions, and in the framework of different macro rheological models 
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(Newtonian, Power Law, Casson, Quemada), our first objective is to characterize 
spatially variable hematocrits and describe quasi-stationary velocity fields and their 
associated energy dissipation with approximate analytical methods. 

Our second objective is the digital simulation of stationary and pulsated flows, with 
axial accumulation of blood cells, applying finite elements methods. 

Our third objective is to briefly describe the effects of the hemodynamics in the 
parameters of fluid-to-wall mass transport processes, applying a_ simplified 
mathematical model based on diffusive boundary layer theory. 


(C) Materials and methods 


Blood flow in in straight segments of small arteries is described here by one phase 
incompressible and isotropic theory of viscous suspensions (Slattery 1972). At each 


point 7 and instant+, the strain rate tensor D(t,7) and the stress tensor G(t,7) are given 
in terms of the fluid velocity field a(t,7), pressure P(t,7) and the deviatoric stress tensor 


(extra-stress tensor) T,,(t,7): 


b=+-(vi+vi") (1) et a a a (2) 


As usual, V represents the gradient operator, ( y the transposition operator and T the 


unit tensor. Scalar shear stress and strain rate are defined thus (Slattery 1972): 


pe [i -tracelé, ##,) 3) j=\2-maceled) 


A general constitutive equation for blood, considered as an isotropic incompressible 
fluid with apparent viscosity iH * 7) and local hematocrit H, can be written: 


~ 


Zp, =2-u(H,,7)-D (5) 


From the scalar shear stress and strain rate, and equation (5), the scalar constitutive 
equations for an isotropic fluid are obtained: 


t= W(H,,7)-7 (6) y= G(t,H,)-t (7) 


If the scalar shear stress and strain rate are known in each point and instant, and if the 
apparent viscosity is known as a function of the local hematocrit and the scalar shear 
rate, from (6) it is possible to determine the local hematocrit field H les r ). 


The one phase fluid dynamics equations for an isotropic and incompressible fluid are 
(being p blood density, a volume weighted combination of plasma and cells densities): 


p-(Zrieva|--vP 0-4, + pF (8) V eu =0 (9) 
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Neglecting the term due to external forces p-F , from equations (1) and (5) it follows: 


p:(Esaeva|=-ve+v-(u(i,7)-fva+vi") (10) 

The spatial and temporal evolution of the local hematocrit field is given by the 
H |. 3 

continuity equation: a +ueVH,=-VeJ (11) 


Here the vector field J represents the flux density of blood cells (mainly red cells) 
relative to the bulk suspension velocity uw. From a non-equilibrium thermodynamics 
point of view (Lebon et al. 2010), it can be considered as an internal variable added to 
the external blood variables p ,H,,u and T (the assumed uniform temperature of blood). 


This one phase continuum approach allows a macroscopic description of red cell 
migration through the space-time distribution of a macroscopic field of hematocrit. 
However, it can’t cope with the cell free marginal plasma layer, adjacent to the inner 
wall surface, which is a microscopic phenomenon relative to scale of the continuum 
model. So the distribution of the local hematocrit at the wall of the conduit, given by the 
mathematical model, should be interpreted as the hematocrit at the inner boundary of 
the cell free plasma layer. 


For concentrated suspensions, like blood, the migration of particles may be suitably 
described by hydrodynamic diffusion theory (Davis, 1996), employing Phillips model 
(Phillips et al. 1992). As shown in van Dinter et al. (van Dinter et al. 2013), the formula 
for particle’s migration current density can be re-written as the product of a positive 
scalar conductance field known as the hydrodynamic diffusion coefficient D,, by 


minus the gradient V of another scalar field I(t, 7), known as the migration potential: 


J =F up =—Dyy *VU (12) 
Here D,,, is proportional to the square of the local particle concentration (the local 


hematocrit H if the fluid is blood), the square of a characteristic particle sizea, and a 


positive measure of the local shear rate in the fluidy: D,,) =K,-H e ay (13) 


The migration potential field is the logarithm of the product of the local particle 
concentration, the local shear rate measure, and a power of the local fluid viscosity yw: 


I =In(H, -7-u*) (14) 


K 
Here A = —is the quotient of two positive non dimensional parameters that take into 
account friction and electrostatic forces during the interaction of charged particles and 
are of the order of unity. Streamline curvatures can be incorporated in the migration 
potential for analytical purposes (Shauly et al. 2000). 
H, . 
From (11) and (12) we obtain: i eVH, =Ve(D,,-VU) (15) 
t 
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Working with equation (14), the following formula for the gradient of the migration 


potential is derived: VII= re pecan ne -VH, + ee een ey ‘Vy (16) 
H 0H y oy 


From (12) and (16) the flux density of suspended particles due to hydrodynamic 
diffusion can be decomposed in the sum of a convective migration term and a diffusive 
term: J =Jip =H,-ii, -D-VH, (17) 


The migration velocity uw, and the diffusion coefficient D are given by the formulae: 


Foe A ee ea oe (18) 
olny 

DIR er (eg (19) 
OlnH, 


For dilute enough suspensions hydrodynamic diffusion ceases to be the main particle 
migration mechanism. However, although now with a different meaning of the 


migration velocity uw, and the diffusion coefficient D, the particle flux density J can 


be written also as J =H, -ii, —-D-VH, (Scott 2005). 


From equations (15) and (17) we derive: ” 


H ee 
< '=-Ve((i+i, )-H,)+Ve(D-VH, ) (20) 
t 
Near the wall, the field of scalar shear strain rate, given by equation (4), influence fluid 
to wall mass transport in the diffusive boundary layer adjacent to the endothelium 
(Ross-Ethier 2002). This has particular importance when the limit step to mass transport 
and reaction is in this diffusive layer, as seems to be the case of oxygen, ATP, 


bradykinin, angiotensin, and other effectors degraded or chemically modified at the 
surface of the endothelium (Caro 2001, Tarbell and Qiu 2006). When the chemical 
species are carried by blood cells, as is the case of oxygen, the lowering of the local 
hematocrit near the wall, due to axial migration of blood cells, could affect the 
concentrations of the dissolved substance and decrease fluid to wall mass transport 
(Biasetti et al. 2014). 


(C.1) A framework for a simplified analytical approach 


To ease the application of analytical methods, a simplified model of axisymmetric 
flow in a straight segment of artery is considered under quasi-stationary conditions. 
Cylindrical coordinates r,@,z will be used. For transient conditions and a cylindrical 


conduit with rigid walls, the fluid velocity is parallel to the axis z , equation (5) and (10) 


* Hematocrit kinetics can be described as a continuous random process if the local hematocrit 
Hi, (t, 7) is interpreted as a probability density of finding a red blood cell in each position in a 


given time instant. By its appearance, equation (20) seems to be closely related with a Fokker- 
Planck equation (Gardiner 1986). However, its coefficients are functions of the probability 
density H, (t, 7). 
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reduce to the following (being L the tube length) and|..| the absolute value: 


clord=af eon ]er]) or) 
pce 2A) Ne eee) (22) 


Ot L r or 


The time varying pressure difference between the inlet and the outlet can be expressed 
as its maximum AP,, multiplied by a suitable function of time: AP(t) = AP, f (t) (23) 


Now let us consider quasi-steady conditions to calculate the pumping power 
E=AP-Qneeded to sustain a given flowQ. Under quasi-steady state conditions, 


equation (22) simplifies to: AP() z : : g (r . t(t,r)) (24) 
L r or 


If R is tube radius, let us introduce the dimensionless radial distance: y= 7 (25) 


From equation (24), we obtain the (well known) rheology independent formula for 


shear stress (Caro et al. 2012, Chandran et al. 2012): c(t, r) = Rar) -y (26) 


To calculate flow as a function of AP we begin with the formula for shear rate as a 


function of shear stress and hematocrit: = g(z,H,)-c =h(z,H,) 


éu(t,r) 1 du(t, y) du(t, y) [4 - AP(t) 
N y= =o —~-— = R-h| ———“ - y, H(t, 
ow AG r) Be 2 By and ay ae y At y) 


R 1 
Blood flow is given by Q(t) = [2x rT: u(t, r): dr = 27 - R* ‘Ly : u(t, y): dy. Taking into 
0 0 


account the non-slip boundary condition at the wall u (t, y= 1)= O and integrating by 


1 
parts it follows: Q(t)=7-R?-| 4( aro ‘y, 11(.7)| -y?-dy (27) 
, ; 
: 3 AP 
From equation (27), we can derive: ae =f (Q,R) (28) 


Then we can estimate pumping power to sustain a given blood flow in the conduit: 
ExQ-AP=Q-L- f(Q,R) (29) 


To be able to use equation (29) it is necessary to determine 7 =h(c,H,) by means of a 
rheological model that considers the effect of a varying hematocrit. For the present 
purpose, the validity of the known relations between the parameters of popular models 
of blood macro rheology (Newtonian, Power Law, Casson, Quemada) and the global 
hematocrit H,,is extended to a relation between the local parameters and the local 
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hematocrit H,. This assumption seems to be a reasonable approximation when the 


derivation of the relation between the viscosity and global concentrations (Lapasin et al. 
1996, Starov et al. 2002, Scott 2005) is re-examined. 


The constitutive equations for Newtonian models, with a local viscosity dependent of 
the local hematocrit (that could be called a generalized Newtonian models), are: 


(t.r)= gy(H,(tr))- Ht.) (30) 


Here we assume, being yw, plasma viscosity, H,,, a maximum packing fraction (often 


t,M 


taken as 1) and k a dimensionless parameter (often taken as 2): 


dy (H, > 2 G1) 


H k 
= t 
| Ay 


The non-Newtonian Potential Law model is given by t= jp,(H,,7)-7 with an 


apparent viscosity Mp, (H,.7) = K (HL, ). y""' . For blood, to describe shear thinning, 


the positive parameter n is less than one. It can be derived from a fractal model 
(Lapasin et al. 1996). 


The non-Newtonian Casson model is given by T= u(H,,7): y with an apparent 


2 
viscosity: u-(H,.7)= | 2" 4K (HA, } (32) 
Whenr < 3 the shear strain rate is zero. Whent > tts the shear strain rate is an 
increasing function of the shear stress. From (32) it follows that if y grows without 
bound, the asymptotic viscosity is K.(H j= Ly (H +) (Here we suppose that ly (H ) is 
given by equation (31)). The yield stress z,(H,) is an increasing function of the 


hematocrit. As the hematocrit can vary from place to place and from instant to instant, 
the yield stress is in general a field t : (t, r ). 


The non-Newtonian Quemada’s model has a constitutive relation t= U(H,.7)- y 


with an apparent viscosity (Quemada 1978 (a) and 1978 (b)): 


U(H,,7)= (33) 
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Intrinsic viscosity [u] is defined by [u] (7)= jim = {2 (H,.7) ] and in 
30 HL, 


ul +b. | 


Quemada’s model it is given by: [u] Ga) = (34) 


Now, from the constitutive equation x(t, y)= u(H, (t, y), Wt, y)) Wt, y) of one of the 


rheological models described above, from formula (26) for the rheology independent 
shear stress, and from a suitable analytical formula for the shear strain rate, it is possible 
to derive a formula for the local hematocrit. 


Let us consider this last topic. Axisymmetric flow velocity fields in small arteries and 
arterioles can de described by an empirical formula (Popel and Johnson 2005, Dillon, 
2012): 


u(t,y)=u, -(l-y*):f()  O<y<1 (35) 


Here y = a fe (t) is a periodic function of time, u,, is the maximum velocity on the axis 


and q = 2 giving the degree of blunting of the profile. In this case the flow is pulsatile, 
not oscillatory, so f (t) is always positive, with maximum equal to | and minimum 


greater than a lower bound of 0.3. It is phase-shifted relative to the cardiac cycle, 
depending of the position of the conduit in the arterial tree. 


1 
The pulsatile flow is: Q(t)= 27 - R? ‘fy-ule, y):dy=a-R?-uy, 5 fo (36) 
0 q 


; 1 . 
The shear velocity is given by: Ht, y) = ao ae y) =q: u yt. f(t) (37) 
R oy R 
: R-AP,, 
Formula (26) for the shear stress can be written: r(t, y) = aap yf (t) 


So both, shear rate and shear stress are known in each point of the conduit. Substituting 


the known values of r(t, y)and Wet, y) in the apparent viscosity formula u(H a. 7) = an 
v 


an approximate formula that allows the determination of the hematocrit field is 


u R? - AP, 
obtained PB MLY ge yt) |e (38) 
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Because the assumed shear velocity is independent of the axial coordinatez, a 
field H,(t, y)independent of z is derived from (38). 


An alternative and more rigorous approach to the derivation of hematocrit field employs 
the migration potential (Phillips et al. 1992, Chen et al. 2003). In a quasi-steady state, 
after the suspension flow is well established, the migration potential II must be 


constant. From equation (14), this is equivalent to the constancy of H,-7- u’(H,,7). 


However, y= and from equation (7): a oe So when there is 


Ll g(t.H,) 


equilibrium in the hydrodynamic diffusion process and the shear stress field is known: 
H,-t-g'* (H,, t) = const. (39) 


As said above, in the case of oxygen, both fields near the wall, of shear rate and of 
hematocrit, could influence fluid to wall mass transport in the diffusive layer adjacent to 
the endothelium. In a first approximation, the flux density J, of a certain substance 


through the diffusive boundary layer can be described by this formula (Tarbell and Qiu 
2006): 


J5= Kew’ (C, an Cy) (40) 


Here C;, is the concentration in the lumen just outside the diffusive boundary layer, C, 
is the concentration in the surface of the endothelium, and k, , is a mass transfer 


coefficient. Assuming a well-mixed entrance flow and a developing diffusive boundary 
layer, the following proportionality relation can be approximately applied to estimate 


the mass transfer coefficient (Tarbell and Qiu 2006, Caro et al. 2012): k, y «3/7, 


Taking into account that 7 =h (z, H “Ns the permeability of the boundary layer verifies: 


kpw ©ih(c, .H,,,) (41) 


It is reasonable to expect an increase in the local viscosity when the local hematocrit 
increases. So, for a given shear stress in the wall 7, (which is independent of the 


rheology and as consequence of the hematocrit and its distribution), the mass transfer 
coefficient should be a decreasing function of the local hematocrit near the wall. ° 


3 In the case of oxygen, it is possible to guess that the concentration just outside the diffusive 
boundary layer can be considered as an increasing function C, (H, ~ of the wall hematocrit. 


The hematocrit field is rheology dependent, as will be seen below. Then it is possible to expect 
that oxygen concentration near the wall is also rheology dependent. 
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(C.2) A framework for a simplified numerical approach 


Velocity profiles, pressures, local viscosity fields, shear rates and wall shear stresses are 
obtained by digital simulation, in straight segments of small arteries under pulsated flow 
conditions. Hematocrit fields are introduced in the viscosity of generalized Newtonian 
models to describe the effect of axial accumulation through viscosity fields in the 
equations of fluid dynamics. 


Due to the symmetry of both the spatial domain and the flow, finite differences could 
have been used. However, finite elements (Ames 1992, Gershenfeld 1999) were 
applied, to ease the comparison with more realistic settings in which less symmetric 
spatial domains appear and spatially inhomogeneous discretization grids are more 
convenient. 


Well known commercial software (COMSOL) already used by one of the present 
authors (Suarez-Bagnasco et al. 2013) was applied to the digital simulation of the 
equations. After a conversion to a weak form, the transport equations are solved using a 
version of a (Galerkin based) finite element discretization method. Possible spurious 
oscillations are controlled by numerical stabilization techniques available in the digital 
simulation code, without having to recourse to an impracticable (due to hardware 
limitations) dense mesh (COMSOL 2011). 


The effects of the resultant hemodynamics on the parameters of fluid-to-wall mass 
transport processes are estimated from wall shear rates by means of approximate 
analytical formulae of diffusive boundary layer theory (Tarbell and Qiu 2006). 


(D) Results 
(D.1) Analytical results 


Quasi-steady flow fields 


From equation (21) and (22), introducing an arbitrary viscosity field ult, r) it follows: 


Ou(t,r) _ AP(t) 1 Of. Ou 
Oe Og of, u(t, r) (.n) (42) 


Now, let us make a change of time and space scales and transform equation (42) to 
apply perturbation methods in the same spirit of the book by Lin and Segel (1988). 


: : : : r : : 
Introduce the dimensionless radial distance y = eo to tube’s radius R, the 


: : ; t F ; F 
dimensionless time fT, = a relative to the time scale 7, of the pulsatile pressure, a 
P 
2 
. EL, . . 
= — (beingv, =— the kinematic 
P Vv, P 


py R 


plasma viscosity diffusion time scale T,, = 
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T, 
viscosity of plasma) and the dimensionless parameter ¢ = cae z Then (42) can be 
Pp P 
R’ - AP 
rewritten: é- u(t» y) = (r,) te ye He p»¥) a (c,,y)| (43) 
Os peek y Oy L, Oy 
p R 1 p R : ee 1 
Here ¢= ‘+ -——-f, being (by definition) f, =— the frequency of 


p Pp v Pp T, p Pp Vv Pp Tp 


flow variation (heart frequency in Hz). Now introduce the Womersley number for 


plasma: a,=R- mua (44) It follows that:  ¢= ores, a, (45) 
V, 2m Pp, 


For blood © ~1.024 and v, =1.16x10°"/ Then, for an arteriole with R = 300 sm 
Pp 
ee ae R’ e 
the plasma viscosity diffusion time scale T,, = * 7.76x10™ s 
Vv 
Pp 


0 (cp, y) 


If the dimensionless number ¢ is less than 0.1, the inertial term ¢- 
Tp 


neglected in a first approximation. For the abovementioned arteriole this happens when 
T, 1s greater than 0.795 s (a frequency of the pressure wave less than 1.26 Hz). 


may be 


; ; ee 
Figure 1 shows the relation between frequency f, = ra (in Hz) of the pressure wave 
P 


2 
and the arterial radius R (in m), fore =0.1: ~ =1.26- (eae? (46) 


P 


0 0,002 0,004 0,006 0,008 0,01 0,012 
R (m) 


Figure 1 Estimation of an upper bound of pulse frequency (to be able to neglect inertial 
effects) as a function of the artery radius for radii between 300 um and 1.1 cm. 
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Quasi-steady hematocrit fields 


If particle migration is produced mainly orthogonally to the fluid velocity field, 
considering the assumed symmetry of the flow, in cylindrical coordinates, equation (15) 


simplifies to this one: ea at ye ie Jn) (47) 


Introducing a representative fluid velocity u, (the velocity on the flow axis), tube’s 
radius R , the global tube hematocrit H, ; and a measure of particle sizea, a suitable 


order of magnitude estimation of the hydrodynamic diffusion coefficient 


Diy =K.-H,-a?-7 is: Dz. ~H, wae (48) 


: . : : r . : 
Introduce the dimensionless radial distance y = Re to tube’s radius and the 


‘ ; ; t ‘ : ‘ 
dimensionless time Tt, = a relative to the time scale T,of the pulsatile pressure. 
P 
2 


Define the hydrodynamic diffusion time scale T,,, = and the dimensionless 


@,* 
3 
number: = ED 7 5 (49) 
T, 12 Oe ‘a tl, 
Then equation (43) can be re-written: 
H 
Ae eny) _ 1 g y Pw. ees) (50) 
OTp y Oy Dg» Oy 
6H, (cp, y) 


Iféis small enough, and if remains suitably bounded in the time scale7,, 


Tp 


equation (45) can be approximated by its steady-state version: 


1 a [> Pe é nls,.»)|=0 (51) 
y oy 


Then the function oie. JO Css must be (practically) independent of the 
y D.. od pry p » p 
@,» OY 


dimensionless radial distance. As it is reasonable to suppose that 5 Ulee. y) remains 
y 


bounded when this radial distance approaches zero, it follows that <n, y) x0. So 
y 


in a quasi-steady state the migration potential must be (practically) independent of y . 
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Axial accumulation of formed elements and blood pumping power 


Now let us estimate the pumping power, given in general by formula (29), for different 
rheological models (Newtonian, Power Law, Casson), considering axial accumulation. 
As already said, the flow in the conduit is restricted and the pressure head is adjusted to 
maintain a given flow. So, the first step in the procedure is to calculate flow according 
to formula (27). The second step is to derive the pressure gradient as a function of flow. 
The third and last step is to calculate the energy dissipation according to (29). 


-Generalized Newtonian model 


In this case y=hy (c,H,) = maa T Ih (H, (t, y)) = My (t, y), formula (27) 
Hy 


a-R*- AP(t) eee AP(t) 5 
gives: Q(t) = ).f y: 3) dy (52) 
0 My ( f, 
From formula (52) the pressure gradient across the vessel follows: 
AP 2-Q 1 
= ,R)= ; Be 
Ms eras = 4 (53) 
| ue dy 
0 Hy (t, y) 
Then, for the generalized Newtonian model the dissipated power is. 
pee, ! OQ? (54) 
a-R r 3 1 
| y?- -dy 
0 Hy (t »y ) 


1 
ree 1 
The dissipated power decreases when F [ v= | ara ee: increases. The 
0 Hy\t.y 


functional F [en | increases when the Newtonian viscosity field decreases towards the 
vessel’s wall. The viscosity field ,(H, ) decreases when the hematocrit H, decreases. 


So, in the framework of a generalized Newtonian model and equal other conditions 
(ceteris paribus), axial accumulation decreases energy dissipation and under isothermal 
conditions decreases entropy production. 


-Power Law Model 


For the Power Law model: vy =hp, (z, H,)= O0<n<l 


wy 
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If K,,(H, (t, y)) = Kot, y) formula (27) gives: 


it eo 
: [K p(t, y)] 
AP 2-Q" 1 
Then: — = f(0,R)=—=_ - (56) 
fs o e 7 ; a4 1 

bee 14 
: [Kpi(t, y)] 

From (29) the dissipated power for a given flow is: 

Pps 7 _.9""(t) (57) 


x" .Re7 r ,; 
Jy *: ——— 4 
: (Kole) 


2 
n 


1 
The dissipated power decreases when F [Hor |= =[ y —————- -dy increases. The 
0 


[k PL (t, y. a 
functional Fi,, | increases when K,,(H, ) decreases towards the vessel’s wall. But the 
field K,, (H , ) decreases when the hematocrit decreases. So, also in the framework of a 


Power Law model, and under equal other conditions, axial accumulation decreases 
energy dissipation and under isothermal conditions decreases entropy production. 


-Casson model 


In this case the shear rate vanishes when0 <7 <z7,, and whent 27, the shear rate is 


ir (H,) 2 
given byy =he. (c,H,)= Wve ry (H, ] 

Hy (H,) 
that appears in the denominator is an asymptotic, hematocrit dependent Newtonian 
viscosity. The yield shear stress is assumed to be given by the empirical formula 
(Merrill 1969, Quemada 1978 (b)) where2 < p<3andH,, =0.17: 


t= 


. The viscosity field (42, (H,(t, y)) = y(t, y)) 


r,(H,)=7,(t,y)*c-(H, -H,,.)? (58) 


“Under in vivo conditions always H, > H, , 
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Wee) )- Jz, ( (i 


Equation (27) gives the flow: Q(t) =7-R?>- fr -dy (59) 
4 Hy(t, y) 
Introduce the auxiliary variable: x= — (60) 
From (59) and (60) it follows: 
5 
1 3 y?. | (t, ) 
ar) (i ay) 2: (2 iad -dy af fo v(t y) ) (61) 
» y(t, y) > My(t.y) ae ss, 


Now, define a functional of the Newtonian asymptotic viscosity and the yield stress 
5 


fields: 


(62) 


1 3 
(5-4 
0 


Then, define a functional of the Newtonian asymptotic viscosity and the yield stress 
fields that is also a function of the blood flow and vessel’s radius”: 


al) fy. ty (ty) i 


m-R J Uy(t.y) 


¢ [uy.7,;(2.R)]= (63) 


The solution with physical meaning of equation (61) can be written thus: 


x=b+,/b’ +c (64) 


Taking into account (60), the pressure gradient to maintain a given flow is: 


= f(Q,R)= z 2 64 r+) (65) 


AP 


> From of or = and the behavior of the integrands there it 


aera ae -dy : the functional defined by equation (63) takes always 
Eto) y y eq 


positive values. 
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Then the dissipated power for the Casson model in a quasi-steady state is: 
. 2 iz L 2 
BOA Ne +c}-Q (66) 


If the hematocrit H, decreases when the dimensionless radial distance increases, both 


the asymptotic viscosity ,(H,) and the yield stress Ty (H,) (formula (58)) decrease 


3 


y 
y(t. 


y “Ir (ty) = ye . T, t,y 


My(t, y) ~ My(t.y) x 


from the vessel axis towards the vessel wall. Then the function 


increases 


Nin 


and 


towards the wall. However, both 


3 3 
2 Ey (y) eee a (y) can’t increase as fast as ——~ when the radial 
pity) alts) 9 My(t, y) 


distance increases, or even they can decrease or vanish. As consequence the functional 


b+.,/b* +c decreases if the hematocrit decreases from the vessel axis towards the wall: 


then, formula (66) shows that also in the framework of Casson model energy dissipation 
should decreases with axial accumulation. 


Hematocrit fields from q-adic velocity profiles: a derivation of model m 


The following mathematical formula for a spatially variable hematocrit (model m ) was 
proposed for axisymmetric velocity patterns in straight segments of small arteries under 
quasi-stationary conditions, as a_ function of the dimensionless _ radial 


, H A, Ww 
distance y = ” (Sudrez-Antola, 1980): Ly) =1-/1 sad (a (67) 
R A, y 1M 
R 
23 Ht [r-A() dr i 
From the tube global hematocrit H,, = re =2: } y: H,(y) -dy a 
1: : 
F P (H, -_ H, ) 
formula for the parameter m is derived: m=2- “8 (68) 
A ay a ca 
If mis known, wall hematocrit can be estimated: H,,, =H, — 3 : (H,,, “a H,,) (69) 


Figure 2 shows the parameter m as a function of the tube global hematocrit ° for the 


® The discharge hematocrit is obtained collecting a volume of blood circulating through a vein 
in a suitable recipient. If the velocity profile u(r) is steady and axisymmetric, and if Q is the 


R 
2-2-[r-u(r)-H,(r)-dr 
volumetric flow in the vein, the discharge hematocrit is: H,,, (disch) = 0 


Q 
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caseH,,, =OandH,, =1. 


0,1 0,3 0,5 0,7 Ht,g 0,9 


Figure 2: Relation between m and global hematocrit (for H,, ranging from 10 to 95%) 


Figure 3 shows the function H Ay) with = for five values of the global tube 


hematocrit for the particular case H, ,, =QandH,,, =1. 


—Htg + 20% 
0,300 —Ht ge 35% 
—Htg + 45% 
— iit 55% 


0,100 —=HtEE 80% 


0 0,2 0,4 0,6 0,8 1 


Figure 3: Local hematocrit as a function of distance from vessel axis for different global 
hematocrits Hi. 


The tube hematocrit can be measured stopping the flow and emptying the vessel blood contain 


R 
2--fr-H,(r)-dr 
in a suitable recipient and can be calculated by the formula: H,,, (tube) = 0 


mR? 
(Baskurt et al. 2007; Nichols et al. 2011; Chandran et al. 2012) 
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Now, let us explore the possibility of deriving this modelm from known quasi-steady 
velocity profiles. 


For quasi-steady axisymmetric flow velocity fields that can de described by the 
empirical formula (35) (u(t, y)=uy (1-y*): f (t)), the hematocrit can be obtained 


from equation (38) because in this case both shear rates field Wt, y) =q: ao yas (t) 


R-AP, 


ap -y- f(t) are known. 


and shear stresses field r(r, y) = 


The rheology independent field of shear stresses can be related with the pulsatile wall 


Pu é * = . 
a EOD (70) 


shear stress 7,,(t) by: t(t,y)* 


The rheology dependent field of shear rates can be related with the pulsatile wall shear 


rate 7,,(t) by: y(ty)=4- yt fQ=7,.0-9" (71) 


To continue, it is necessary to specify a formula for WH,.7), that is a rheological 
model, in formula (38). Let us introduce the viscosity of blood adjacent to the vessel 


R?-AP 
wall: dice fw = Ml (72) 
%y 2°q- Letty 


Taking into account (70), (71) and (72), formula (38) can be rewritten thus: 
HH, (ty). A) y= My (73) 


Quasi-steady hematocrit field when the apparent viscosity is given by the 
generalized Newtonian model 


If the viscosity is given by formula (31) ( Ly (H ————\—_ ), from (73) it follows 
i ay 

. A, Ww m 
model m formula for the local hematocrit: H, (y) = Ay (. [- }> | (74) 


= —2 A, _| Hp 
oor ia (q2=2) (75) t as |_(2| (76) 


From (74) and the formula that relates the local with the tube hematocrit 


1 
H,, =2-[y-H,(y)-dy , the wall hematocrit H,,, = H,(1) results as function of 
0 


maximum packing hematocrit H,,,, the global (tube) hematocrit, and the parameter 


m (determined in terms of the assumed known parameters g and k from (75)): 
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A, =»: ~7.(H,,,-H,,) (77) 


H H : 
From (77): 1-—" =| 1-—*£ (2 From (75); 14 =144°% (g>2) 
t,M Ay 2 2 pi 


From the formula that relates the viscosity with the hematocrit in a generalized 


= ; and dy (H,. )= wy. = EE 


H, H,,) 
i= [ea eee 
Ay Ay 


Then, from these last four equations we obtain a formula that gives w,,as function of the 


Newtonian model sy, (H ie ) =U, ® 


parameters ww, kand H,,, of the relation between the viscosity and the hematocrit, the 


blunting parameter gq of blood flow profile and (through 4, ,) the global hematocrit 
k 


Hs i, © a) ‘Hy (78) 
t,g w = 74 
1+ ine 


Quasi-steady hematocrit field when the apparent viscosity is given by 
Quemada‘s model 


If the viscosity is given by formulae (33) and (34), 4, (H i. y)= 


eb +L. fz 


[ul(y) = ——_—_—__——. , from equations (73) and (71) the following formula for the 


local hematocrit is obtained: 


{vs Flt | 1 
k 
H(t, y) = = ih ta “y" (79) 


[U4 Ual.- | ye | 


The parameter m is also given by formula (75) as in the generalized Newtonian model. 
However, instead of a constant maximum packing fraction H,,,, a field H, , (t, y) of 


local maximum packing fractions results: 


(80) 


cca alas | 
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. : : ona ae : 
When the maximum of the dimensionless shear rate tolt) is small enough, or when its 


Ye 
minimum is big enough, formula (80) gives a constant maximum packing fraction, 


one and gee: respectively. 


leh Lal. 
Asymptotic steady hematocrit field when the apparent viscosity is given by 
Casson’s model. 


Let us consider strictly steady conditions and describe blood rheology with Casson 
model. When the shear stress 7,,- y at the dimensionless radial distance y is greater 


than the yield stress 7, the shear velocity there is different from zero. The yield stress is 


a local value that increases when the local hematocrit increases according to equation 
(58). Here we suppose that H (y) > H,,,for every dimensionless radial distance. When 


there is axial accumulation 1, (H,)= Ty (y) decreases from a maximum value 7, (H,,) 
in the vessel axis towards a minimum value at the wallr, (H,,,). Local shear stress 
t(y)= T° y imcreases from zero in vessel axis towards 7, at the vessel wall. 
Ifz,, >t, (H,,,). the decreasing function 7, (y) and the increasing function t(y) have 
one intersection point y. given by the unique solution of the equation t(y) =, (y). AS 
we don't know yet the field H (y)we can’t specify y., but we know that there is one, 


and we know that when y, < y, theny > 0. In this case equation (73) gives: 


2 
T (A, m 
alt de-s el [OO «ast =p,-y" (81) 


However, if y<y. then y=0 (we have piston flow in the inner cylinderO < y < y, ). 


Hp 


Taking into account that /u, (H,) = ;» equation (81) can be recast as 
H 2 
A, y 
q-2 
P 2 
—— ‘y 
7(H,) 1 Vu, 
follows: ‘at pal y2y, (82) 
Mw 2 Hw \2 
sy, 
; (ane oer ; 
Now, let us define ¢ = ./~~—*~— so (82) can be written thus: 
qT, 
Pisogiat 
TAA 1 
ra eA ee oe (84) 
eg) 5 2 
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TAH uw). : 
———— is small enough, a perturbation 


w 


When the dimensionless parameter ¢ = 


analysis gives the following successive approximations: 


k 
Zero order (€ =0): H, (y)* ay: 1-[ 4 “y's (85) 
oe 


k 
q-2 H, ii, x 
Here m= i and =|1 So, under steady conditions model m formula 


/,, t,M 
is obtained as the zero-order approximation. 
TA). bots ; . 
First order (€ = ,| ——— is retained in equation (84) and 1, (H, ) is substituted in 
ts : : 
place of 7, (H, )): 
q 
P 2 
pale y 
Ty H, oly 1 L,, 7 A é a 
() + 71 From this approximate equation a first 
- is 
Ww 2 


a oy Ty H,,(y) 1 
l-e- ; 
r., Hy) ‘2 
Now, make the further approximation: 
H 
d wlt+e- rs tal) : : 
T, A, y 1 *y - y? 
é . = . 
\ t,(H, 1) y2 
Then this formula for the hematocrit field is finally derived: 
1 1 
ML, \k Hye [ty(Hioly)) m4 
H,,(y)* Hy *| 1 a (2 Ooi ee (86) 
af ) - [ee LL, T,, ey 


When m= t= > ; or equivalently when g>2 + , the decreasing function of y 


given by formula (86) is well behaved in the whole interval of radial distances. 
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Hematocrit fields from the migration potential and from Scott’s micro- 
rheological model 
Hematocrit from migration potential 


When there is equilibrium in the hydrodynamic diffusion process and the shear stress 
field is known, equation (39) gives the corresponding condition of equal migration 


potential everywhere in terms of the function f= (c,H = a introduced in 
r w(H,,7) 
equation (7). From (39) and (70) the following formula is obtained for the hematocrit 
A-l 
H,,,(t) [ae (87) 
A, (t, y) o(H,,, (¢), Ty, : f(t) 


The parameter 4 is greater than | and known. The relation between tube and local 


field under quasi-steady conditions: y= 


1 
hematocrit (H,, = 2-[y -H, (t, y)-dy< H,,) allows the derivation of the wall 
0 


hematocrit for a slow enough pulsatile flow. Wall shear stress 7, is known and 
rheology independent. The positive function f (t) is assumed known. The size of the 


suspended particles doesn’t appear in equation (87). So, in principle, the quasi-steady 
hematocrit field predicted by the hydrodynamic diffusion theory is independent of 
particle’s size. 


k 
-For a generalized Newtonian model, (ett,)=E-| 1 at . In this case (87) 


Hy 1M 
0)" 
gives: y= Bow Ao (88) 
A, (y) t= A, w 
t,M 


1 H, ‘ 
pei) (89) 
k- t Er | 
Here, by definition: 13 Ba (r)= = (90) 
lb 


In (90) 7, is a constant critical shear stress and t= 7(t, y)=7,-y- f(t) 


25 


Research Report September 2014 


H, (i, y) (A-1)-k 


Hult) | Aywley-y- FO) 
. _ t,w , t,M Ty : y 
Then, from (87) and (89): y= H, (. > He, i) (91) 
A, M qT, ; f(t) 
2 
H H 
-For Casson model, whenever rs 2 a) y At, H,)= a 1 rl “| and 
T Ly (H,) T 
ict ply ay) 
y At.H,) is zero if“ > 1. From (87) and Ty (H(t, y))<t,,-y- f(t), it follows: 
c 
; Hi (A-1}k r, ( A, (A-1)-2 
A, Ay ‘ ty f(t) 
: A, _ Ay 1 CAH, a 
Ay T 


Under strictly steady conditions, we have seen above (in the discussion immediately 
before equation (81)) that, (H (y))= T,,:y has only one solution, the critical radial 


distance y,, such that there is piston flow in the inner cylinder O< y< y. and in the 


(A-l}k (A-1)2 
= A, ia Ty t 
H H T° 
annular region outside it: y=": = ries (93) 
H,(y) {= H,,, i 7, (H, 
Ay © 


For blood k is always very near 2. As first approximation it is possible to take 2 = *. sO 


ifk =2 then (A = 1). k =1. This allows the derivation of explicit analytical formulae for 


the hematocrit field. For the generalized Newtonian model and the modified 
Quemada’s model the problem reduces to an algebraic equation. For Casson model, if 


7 \H 
é= j——"- is small, it is possible to obtain analytical formulae by a method of 
Tq, 


successive approximations, like the one already applied for the derivation of an 
asymptotic steady hematocrit field, in the framework of Casson’s model and q-adic 
velocity profiles. The corresponding calculation for a constant migration potential will 
not be done here. 


-If Ci - 1). k =1, from (88), for the generalized Newtonian model the hematocrit field is 


given by the following formula, already given in Phillips et al. 1992: 


H 
aM (94) 


H,(y)=—— 
1+|—" -1}. 
‘leah 


The wall hematocrit H,,, is determined from (94) and the tube hematocrit formula 


26 


Research Report September 2014 


1 
A, , =2-[y-H,(y)-dy (95) 
0 


The following implicit equation is obtained: 


H 
: iw 41 : -log,| —* (96) 
tg Se 
Ay Aly H,.,, 
pea Sic ey | pee eee | : 
H 


If (A-1)-k =1, for the modified Quemada’s model, from (90) and (91), the field of 
Ae (r) 


Here H, y (z) is given by equation (90) and the shear stress isz = c(t, y)=7,,-y- f(t). 


hematocrit 1s: H(t, y)= (97) 


The wall hematocrit H y(t) is determined from (96) and the tube hematocrit formula 


giving another implicit equation if the global tube hematocrit is assumed known. 


Hematocrit from Scott’s model 


Under quasi-steady conditions and with an axisymmetric flow pattern, radial diffusion 


D 
and radial migration flux densities will be equal and opposite: H,-u,, = Ph ~ H, (98) 
y 


In general, both particles migration velocity u,, and diffusion coefficient D will be 


functions of the radial distance. Here we assume axial migration, sou,, = Su 


R-|u,,\(y) 
D(y) 


oil 


m 


Let us introduce a local dimensionless Péclet number: Pe (y) = (99) 


Taking into account that H,(y) must be always less than the maximum packing 
fraction H, ,,, to integrate (98) under steady conditions from a dimensionless radial 


distance y up to the wall, two cases should be considered: 


1 
If H,,, a} } Pe(E )- dg <H,,, the hematocrit field is given by the formula, valid for 
0 


everyO< y<l: Ho) Hy. | Pe) (100) 


1 
If H,,,-exp| | Pe(¢)-d¢ |>H,,,, let us introduce the critical radial distance y. such 
t,w Pp t,M 
0 


that: H,, os} | Peles =H y (101) 


1 
Then H,(y)=H,,, os } rele in the interval y, < y <1 (102 a) 


y 
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In the intervalO< y<y,, H,(y)=H,y (102 b) 
In any case, the local hematocrit field H (y) is related with the tube hematocrit by the 
1 
(often quoted) formula: Hig 2: | y-H,(y)-dy<H,y (103) 
0 
To simplify let us assume u,, (y) =-u,,-yand D=D, constant. Introduce the 
dimensionless wall Péclet number: Pe,, = Pe (1) = Ritty (104) 
0 
Then the steady hematocrit field is given, from (100) and (104), by the formula: 
le 
H,(y) = H,,, oxy] . (iz a) = H,(0)-exp| ~2° . | (105) 
; Tigeke ; : 
Then, from (103) it follows H (0) = = P and from this last equation and 
e 
27) L=exp|=—— 
[ : | 2 ) 
Ee oPe. P 
from (105): H(y)= zs 0] . ; | (106) 
2 fi eS exp|- ew } 


If the maximum packing fraction is attained at a positive critical radial distance y,, this 
parameter verifies the equation (derived from (101)): 


Hey : exp] ( = x) = Ay (107) 
Then: H,(y)=H,,, x z(t a) y.<y<1 (108 a) 
H(y)=H,y OS ySy, (108 b) 


> 


gives a measure of the magnitude 


(2) (B-x.) 


of the migration flux relative to the diffusive flux: in principle Pe could increase from 
near zero, for a uniform distribution of cells in a cross section of the vessel 
(corresponding to a diffusion-dominated regime) towards values an order of magnitude 
greater than 1 for an hematocrit highly concentrated in the neighborhood of the vessel 
axis (corresponding to a full migration-dominated regime). 

Scott obtained an approximate analytical formula for the migration velocity which 
relates said velocity with particle size (through a particle radius a ), vessel radius R, a 
measure of suspending fluid viscosity (through plasma viscosity “,), a measure of 


cell’s rigidity (through a surface tension @ ), a measure of the maximum flow velocity 


V 2 3; 
V, and the dimensionless radial distance: u,, = ae ee (2 (4 (Scott 


70 a R 
2005). Then the migration velocity adjacent to the vessel wall is: 
l 1 : V 2 3 
pee a ee (109) 
70 a R 
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As consequence, wall Péclet number increases when the particle diameter increases. So, 
to the extent that this approach can be applied to a suspension of particles with different 
sizes, the smaller particles will be mode uniformly distributed relative to the larger 
ones: in the flowing blood case, platelets will accumulate less towards the axis than red 
blood cells. This is a significant difference with the prediction of the hydrodynamic 
diffusion theory. According to this last theory steady concentration profiles should be 
independent of particle sizes. 


(D.2) Numerical results 


A Navier-Stokes equation (incompressible fluid) with a steady viscosity field (H,(7)) 
(a particular case of equation (10) with a generalized Newtonian model) was solved 
numerically for two rigid cylinders, one segment of an _ arteriole of 
radius R=0.3mmand lengthL=6mm and the other, a segment of an artery with 


R=6mm andL=40mm. A fluid density 105048/ 3 was used in all the numerical 


calculations. 
Due to the symmetry restrictions, numerical solution reduces to the digital simulation of 
fluid flow according to the equation (42) for a steady viscosity field: 


du(t,r) AP(t) 1 0 ou 
omer aia ee u(H,(r)) alter) (110) 
Boundary conditions: u(t,R)=0 , Peak (t)=0, a (t) Several input pulse shapes 


Vie 


input 


shape P,. (t) (Figure 4) obtained experimentally in a hemodynamic work bench (Suarez- 


(t) were introduced, both pulsatile and oscillatory. Two pulsatile pressures: a wave 


Bagnasco et al. 2013) and the sum P,(t) (Figure 5) of an average pressure with a 
harmonic oscillating pressure. Two oscillatory pressures: harmonic oscillations of 
approximately 1 Hz and 10 Hz, P(t) (Figure 6) and P,,(t) (Figure 7) respectively. 


Global: 


440 | i | /\ i i i 4 
420 4 
400 | 4 
380 + \ 4 


Figure 4: Experimental pattern of pulsatile pressure, with systolic value 440 Pa and 
diastolic value 160 Pa, obtained in vitro in a work bench. 
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Analytical formula used to calculate P,. (t ) : 


4*(68.47+29.85*cos(t[ 1/s]*7.666)+1.167*sin(t[ 1/s]*7.666)+5.808*cos(2*t[ 1/s]*7.666)- 
11.41*sin(2*t[1/s]*7.666)+4.27*cos(3*tf 1/s]*7.666)- 

2.542*sin(3*tf 1/s]*7.666)+0.2223*cos(4*t[1/s]}*7.666)+1.579*sin(4*t[1/s]*7.666)+0.2855*cos(5*t[ I/s ]* 
7.666)-0.2226*sin(5*t[1/s]*7.666)-0.3875*cos(6*t{ 1/s]*7.666)-0.0622*sin(6*t[ 1/s]*7.666)) 


Global: 


Pressure (Pa) 


1 
Time (s) 


Figure 5: Pulsatile pressure with systolic value 440 Pa, diastolic value 161 Pa, mean value 
300 Pa and 1 Hz (approximate) frequency. 


Analytical formula used to calculate P, (t ) : 
140*cos(t[1/s]*7.666)+300 


Global: 


Pressure (Pa) 
° 


1 
Time (s) 


Figure 6: Oscillatory pressure of (approximate) frequency 1 Hz, with systolic value 140 Pa 
and diastolic value -140 Pa. 
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Analytical formula used to calculate P, (t ): 
140*cos(t[1/s]*7.666) 


Global: 


fh fA hh = 


Pressure (Pa) 
i=} 
+ 


_ VV 


01 0.2 0.3 0.4 
Time (s) 


of 


Figure 7: Oscillatory pressure of (approximate) frequency 10 Hz, with systolic value 140 
Pa and diastolic value -140 Pa. 


Analytical formula used to calculate P, (t ): 
(140*cos(t[1/s]*76.66)) 


Figure 8 shows one finite element mesh generated by the digital simulation code. 


Figure 8: A finite element mesh for the arteriole. 
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A time interval comprised between 0 a 2s was used for digital simulation purposes, with 
increments of 0.05s. 


Two cases of the generalized Newtonian model were implemented, both with H,,, =1 
and 44, =0.0012 Pa-s. 


Case 1: (= GB (111) 
Case | equation for the relation between viscosity and hematocrit was not considered 
previously in the present work. Here ( = 0.2 is selected for digital simulation purposes. 
In Case 1 a simplified model m was employed: model m(s) that follows from 
equations (67) and (68), putting H,,,=OandH,,, =1: 
m(s A, . 
H,(y)=(1-y""”) a 
Three values of the global hematocrit H, , were used: 35 %, 45 % and 55 %. The 


corresponding values of m(s) are 1.08, 1.64 and 2.44 respectively. 


As consequence of the chosen relations between viscosity and hematocrit, and between 
hematocrit and radial distance, the implemented viscosity field is: 


Hy m(s\(¥) = (112) 


pao) 
Analytical formula used as input by the software to calculate the viscosity field given by equation (112): 


0.0012[Pa*s]/(1-0.2*(1-(sqrt(x*x+y*y)/Rtube)*m)) 


0.0015 F t t t t = 


0,0014 + 


Dynamic viscosity (Pa*s) 


0,0013 + 


0.0012 


1 1 1 1 L 1 n 
10] 25 30 35 
Arc length 


Figure 9: Viscosity fields along an arteriole diameter for different values of ms) 
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Case 2: an(Ht,)=7 a (113) 


1-H,) 
Case 2 equation is obtained from equation (31) with H,,, =1 andk = 2. Equation (106) 


for the hematocrit field was implemented in this case. As consequence of the new 
chosen relations between viscosity and hematocrit, and between hematocrit and radial 
distance, the following viscosity field was implemented: 


Lt, 
Uy sly) = ; (114) 


23) Lex eee : 
2 


Three values of the wall Péclet number were used: 0.2, 1.0 and 2.0 with a global 
hematocrit H,, of 45 %. 


Analytical formula used as input by the software to calculate the viscosity field given by equation (114): 


0.0012[Pa*s]/(1-(A/(1-exp(-A)))*0.45*exp(-A*(sqrt(x*x+y*y)/Rtube)*2))%2 A= (Pew) /2 


0.0029 Ft t + t t t _—— t t t t t > 


0,0028 + 7 x Re 


7k / \ 4 
0.0027 F \ 


0.0026 + / \ 
0.0025 + / ¥ + 
/ 
0.0024 + / 
/ = 
0.0023 + / \ | 


0,0022 + / \ 


Dynamic viscosity (Pa*s) 
—— 
= 


0.0021 + / .. 4 


0,002 + Se a a | 
7] x oS 


0,0019 / \ 4 


0.0018 F f \ 


0.0017- = / \ 


0.0016 | Mt 
0 5 10 15 20 25 30 35 40 45 50 55 x10° 
Arc length 


Pe, 


y 


Figure 10: Viscosity fields along an arteriole diameter for different values of A = 


A whole blood viscosity of 0.0035 Pa-s was employed in numerical calculations of 
blood flow in absence of axial accumulation, for comparison purposes. 
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(D 2.1) Velocity profiles and shear rates for an arteriole 


Input pressure P,. (t), viscosity fields according to Figure 9. 
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Figure 11: Velocity fields along an arteriole diameter for the experimentally determined 
pressure input different values of ms) for systole and diastole. 


— m=1,08, Time=0.8 
—— m=1.08, Time=1.3 | 7 
— m=1.64, Time=0.8 


~~ m=1.64, Time=1.3 
—— m=2.44, Time=0.8 | 7 
—— m=2.44, Time=1.3 


Velocity field, z component (m/s) 


0.005 f i 1 


i 1 
0 5 10 15 20 x10% 
Arc length 


Figure 12: Details of velocity fields near the arteriole wall for systole and diastole. 


Input pressure P, (t), viscosity fields according to Figure 9. 
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—— m=1.08, Time=0.8 
| —— m=1.08, Time=1.3 
—— m=1.64, Time=0.8 
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Figure 13: Velocity fields along an arteriole diameter for different values of m(s) and 
different instants for the considered oscillatory input for systole and diastole. 
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Figure 14: Details of velocity fields near the arteriole wall for systole and diastole. 


Input pressure P(t), viscosity fields according to Figure 9. 
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Velocity field, z component (m/s) 


Figure 15: Velocity fields along an arteriole diameter for different values of m(s) and 
different instants for the considered oscillatory input for different instants. 
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Figure 16: Details of velocity fields near the arteriole wall 


Input pressure P, (7), viscosity fields according to Figure 10. 
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Figure 17: Velocity fields along an arteriole diameter for different values of A = 


different instants for the considered oscillatory input for systole and diastole. 
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Figure 18: Details of velocity fields near the arteriole wall for systole and diastole. 
Input pressure P, (t), viscosity fields according to Figure 10. 
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Figure 19: Velocity fields along an arteriole diameter for different values of A =——* and 


different instants for the considered oscillatory input for different instants 
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Figure 20: Details of velocity fields near the arteriole wall 


Input pressure F,, (t), viscosity fields according to Figure 10. 
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Velocity field, z component (m/s) 
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Figure 21: Velocity fields along an arteriole diameter for different values of A = 
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Figure 22: Details of velocity fields near the arteriole wall 


(D 2.2) Velocity profiles and shear rates for an artery 


Input pressure P, (t), viscosity fields according to Figure 10. 
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Figure 23: Velocity fields along an artery diameter for different values of A = — 
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Figure 22: Details of velocity fields near the artery wall for systole and diastole. 


(E) _ Discussion and Conclusions 


-CFD results show that for an arteriole the quasi-steady solution of fluid dynamics 
equation, with axial accumulation, agrees with the numerical solution when the pulse 
frequency is of the order of | Hz or less. For frequencies of 10 Hz or more inertial 
effects can’t be neglected and the quasi-steady solution is not applicable. 
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-CFD results show that for an artery the quasi-steady solution of fluid dynamics 
equation, with axial accumulation, doesn’t agree with the numerical solution, even 
when the pulse frequency is of the order of 1 Hz. In Figures 15 and 21, notice the 
change in the sign of the curvature of the velocity profiles for time instants in a 
neighborhood of the times in which the input pressures cross by zero. In these instants 
the time derivative of the pressure pulse takes its maximum values, and the 
accelerations are different at different distances from the vessel axis. Blood adjacent to 
the wall, up to the inflection point of the velocity profile, accelerate positively or 
negatively more than the blood in the region adjacent to the flow axis. 


-According to CFD calculations, in the arteriole, when quasi-steady conditions apply, 
systolic and diastolic velocity profiles are the same for the experimental P,. (t) and the 


shifted harmonic P, (t) pressure pulsatile inputs, when pressure maxima and minima are 
the same (440 and 160 Pa respectively). 


-As expected, CFD calculations with model m(s) for the hematocrit field, give the same 
shear rate at the wall for any global hematocrit. The differences appear in this case 
moving away from the wall towards the vessel axis. So, in the framework of model 
ms) the permeability at the boundary layer is not significantly affected by local 
hematocrits. In the framework of model m with a nonzero hematocrit at the wall, the 
permeability is modified when the global hematocrit changes. 


- With £=1 the viscosity field given by equation (112) reduces to suy,,, «(y= a 
yr AY 


A potential law formula like this was introduced by Womersley in equation (110) trying 
to obtain analytical solutions for oscillatory pressure pulses of any frequency. He found 
that for m=2 and m=4 it is possible to obtain closed form solutions in terms of 
known special functions of mathematical physics (Sudrez-Antola 1980). 


-If quasi-steady conditions are not applicable to the haematocrit dynamics in an 
axisymmetric blood flow, the unsteady nonlinear diffusion equation (47) for the 
haematocrit should be solved, suitably coupled with the Navier-Stokes equation, for the 
geometries studied in the present work. However, when there is quasi-equilibrium in the 
hydrodynamic diffusion process, equation (88) shows that the hematocrit field can be 
considered as steady, at least approximately, as was done here for purposes of digital 
simulation. 


-The hematocrit field of model mcan be derived from quasi-steady q-adic velocity 
profiles if a generalized Newtonian macro-rheological model of viscosity is used to 
describe the relation between local viscosity and local hematocrit. Also, model m 
formula for the hematocrit field is obtained from Quemada’s non-Newtonian viscosity 
model, when the shear rates adjacent to the wall are an order of magnitude either 
smaller or greater than Quemada’s critical value. For intermediate values a 
generalization of modelm formula is obtained (equation (79)) giving a time varying 
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hematocrit if the flow is pulsatile. However, model m formula is not obtained from q- 
adic velocity profiles if Casson macro-rheological model is employed, except as an 
asymptotic zero order approximation when the shear stress adjacent to the wall is at 
least an order of magnitude greater than Casson’s yield stress and Casson model reduces 
to a Newtonian model. This last result agrees with the fact that a q-adic profile can’t 
properly describe piston flow in a significant inner cylinder of flowing blood, although 
for gq big enough the blunting of the velocity profile approaches a piston flow. 


-Even without taking into account axial accumulation, the permeability k,,,, of the 


diffusive boundary layer may depend of the rheological model used. For quasi-steady 
axisymmetric tube flow, let us compare the consequences on the permeability of a given 
substance of employing the Newtonian and Casson models. The rheology independent 
wall shear stress 7, will be taken to be the same in both models, and the asymptotic 


viscosity of the Casson model will be taken equal to the viscosity “, of the Newtonian 


We ve)? 


model. Then, the wall shear rate according to the Casson model is 7,, . = —————— , 
Lu 


Tq, 


and according to the Newtonian model 7,, , = 


From these equations it follows: 
N 


, Ty 2 
V%w.c = VV%w.n ~,|—— Then, for the same shear stress7,,, 7,,.i8 always less thany,, \ . 
Hy 


The permeability is a strictly increasing function of 7, (according to a simplified but 
often quoted model, it is proportional to3/7,, ). As consequence, the fluid to wall 


permeability of a given substance, for the samez,,, will be less if the fluid behaves like a 


Casson fluid, relative to a Newtonian fluid whose viscosity is equal to the asymptotic 
viscosity of the Casson fluid. This remark could be used to justify the calculation of the 
approximate fluid to wall permeability (of oxygen and other substances with a transport 
limit step in the boundary layer) from the digital simulation of the effect of blood 
viscosity on hemodynamic factors in a small, bifurcated artery, that appear in the work 
of Kanaris (Kanaris et al. 2012). 


-The function u,, (y)= —u,,: y was derived by Scott applying a theoretical framework 


already developed by Chan and Leal to describe the motion of a single deformable 
particle suspended in a Newtonian fluid. For this derivation the particle was at a certain 
distance from the flow axis, immersed in the parabolic 2-adic velocity field 
corresponding to an already established steady axisymmetric low Reynolds number 
flow in a cylindrical tube. In principle, it can be applied to dilute enough suspensions. 
Nevertheless, Scott argued about the possibility of extending the derivation to 
concentrated suspensions “treating the surrounding particles as a smeared-out effective 
medium” (op. cit. page 46). In such a case, instead of the plasma viscosity, a measure of 
the suspension viscosity would be used and a q-adic velocity profile could be employed 
instead of the 2-adic one. 
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- Hydrodynamic diffusion allows an understanding of a mechanism of axial 
accumulation of formed elements, but this mechanism is dominant only in concentrated 
enough suspensions. In dilute suspensions other migration mechanisms can be 
dominant. Now, hydrodynamic diffusion is closely related with a stochastic approach to 
axial accumulation, like the already classical approach to mathematically model the 
spatial and temporal distribution in population dynamics (Sudrez-Antola 1979, Okubo 
1980). However, as will be discussed elsewhere, this last approach can be formulated to 
embrace concentrate and dilute suspensions. 


To end, summing up: 


-Taking into account spatial variations of viscosity may modify significantly local 
hemodynamics and fluid-to-wall transport processes as well as pumping power needed 
to maintain pulsated flow in a small artery, both with (generalized) Newtonian and Non- 
Newtonian rheology. 


-In the considered global (tube) haematocrit range (35-55%), obtained shear rate 
profiles and transport parameters in fluid adjacent to the wall are affected by the macro- 
rheological model used and by the local haematocrit field according to the micro- 
rheological model proposed. 
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